// 版权归Go作者所有。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// 下面的原始C代码、长注释和常量
// 来自FreeBSD的/usr/src/lib/msun/src/s_expm1.C 
// 并随附此通知。go代码是原始C.
// 
// ================================================================================================
// 版权所有（C）1993由Sun Microsystems，Inc.保留所有权利。
// 
// 由SunPro开发，SunPro是Sun Microsystems，Inc.的一家企业。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// =============================================================================================
// /
// expm1（x）
// 返回exp（x）-1，x减去1的指数。
// 
// 方法
// 1。参数约简：
// 给定x，求r和整数k，使
// 
// x=k*ln2+r，|r |<=0.5*ln2~0.34658 
// 
// 这里将计算一个校正项c，以补偿
// 四舍五入到浮点数时r中的误差。中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网。用
// 区间[0,0.34658]上的特殊有理函数逼近expm1（r）：
// 因为
// r*（exp（r）+1）/（exp（r）-1）=2+r**2/6-r**4/360+。。。
// 我们用
// r*（exp（r）+1）/（exp（r）-1）=2+r**2/6*R1（r*r）
// 也就是说，
// R1（r**2）=6/r*（（exp（r）+1）/（exp（r）-1）-2/r）
// /=1-r**2/60+r**4/2520-r**6/100800+。。。
// 我们在[0,0.347]上使用一个特殊的Reme算法来生成
// 一个r*r中的5次多项式，以近似R1。这个多项式近似的最大误差是有界的。换句话说，
// R1（z）~1.0+Q1*z+Q2*z**2+Q3*z**3+Q4*z**4+Q5*z**5 
// Q1=-1.6666666666666567384E-2，
// Q2=3.968253968137036583E-4，
// Q4=2.5051361420858517002E-7，
// Q5=-6.28435056823826171002E-9；
// （其中z=r*r，Q1到Q5的值如下）
// ，误差范围为
// 5 |-61 
// 1.0+Q1*z+…+Q5*z-R1（z）|<=2 
// | 
// /
// /expm1（r）=exp（r）-然后，通过以下具体方式计算1，以最小化累加舍入误差：
// 2 3 
// r[3-（R1+R1*r/2）]
// expm1（r）=r+-+-*[-------------]
// 22[6-r*（3-R1*r/2）]
// 
// 为了补偿参数约简中的错误，我们使用
// expm1（r+c）=expm1（r）+c+expm1（r）*c 
// /~expm1（r）+c+r*c 
// 因此，将添加c+r*c作为
// expm1（r+c）的更正条款。现在重新排列术语以避免优化
// /搞砸：
// （2）
// （{（r[R1-（3-R1*r/2）]}r）
// /expm1（r+c）~r-（{r*（---[--[-[-]-c）-c}--）英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网英语学习网。缩小比例以获得expm1（x）：
// 从步骤1开始，我们得到了
// expm1（x）=2**k*[expm1（r）+1]-1 
// =或2**k*[expm1（r）+（1-2**-k）]
// 4。实施说明：
// /（A）。为了节省一次乘法，我们将系数Qi 
// 缩放为Qi*2**i，并用（x**2）/2替换z。
// /（B）。为了达到最大精度，我们通过
// （i）如果x<56*ln2，返回-1.0，（如果x！=inf，则提高不精确性）
// /（ii）如果k=0，则返回r-E 
// /（iii）如果k=1，则返回0.5*（r-E）-0.5 
// /（iv）如果k=1如果r<0.25，则返回2*（（r+0.5）-E）
// 否则返回1.0+2.0*（r-E）；如果（k<2 | | k>56）返回2**k（1-（E-r））-1（或exp（x）-1）
// （vi）如果k<20，返回2**k（（1-2**-k）-（E-r）），否则返回2**k（1-（（E+2**-k）-r）
// （vi）特殊情况：
// expm1（-INF）是-1，而
// 对于有限参数，只有expm1（0）=0是精确的。
// 
// 精度：
// 根据误差分析，误差始终小于
// 1 ulp（最后一位的单位）。
// 
// 杂项。信息。
// 对于IEEE双
// 如果x>7.09782712893383973096e+02，则expm1（x）溢出
// 
// 常量：
// 十六进制值是以下
// 常量的预期值。可以使用十进制值，前提是
// 编译器将十进制值转换为二进制值的精度足以产生所示的十六进制值。
// 

// Expm1返回e**x-1，x减去1的基e指数。
// 当x接近零时，它比Exp（x）-1更精确。
// 
// 特殊情况为：
// Expm1（+Inf）=+Inf 
// Expm1（-Inf）=-1 
// Expm1（NaN）=NaN 
// 非常大的值溢出到-1或+Inf。
func Expm1(x float64) float64 {
	if haveArchExpm1 {
		return archExpm1(x)
	}
	return expm1(x)
}

func expm1(x float64) float64 {
	const (
		Othreshold = 7.09782712893383973096e+02 // 0x40862E42EFA39EF 
		Ln2X56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1 
		Ln2HalfX3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73 
		Tiny       = 1.0 / (1 << 54)            // 2**-54=0x3c90000000000000 
		// 与expm1 
		Q1 = -3.33333333333331316428e-02 // 0xBFA111111111110F4 
		Q2 = 1.58730158725481460165e-03  // 0x3F5A01A019FE5585 
		Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7 
		Q4 = 4.00821782732936239552e-06  // 0x3ED0CFCA86E65239 
		Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D 
	)

	// 特殊情况
	switch {
	case IsInf(x, 1) || IsNaN(x):
		return x
	case IsInf(x, -1):
		return -1
	}

	absx := x
	sign := false
	if x < 0 {
		absx = -absx
		sign = true
	}

	// 过滤掉大参数
	if absx >= Ln2X56 { // 如果| x |>=56*ln2 
		if sign {
			return -1 // x<-56*ln2，返回-1 
		}
		if absx >= Othreshold { // 如果| x |>=709.78。。。
			return Inf(1)
		}
	}

	// 参数约简
	var c float64
	var k int
	if absx > Ln2Half { // 如果| x |>0.5*ln2 
		var hi, lo float64
		if absx < Ln2HalfX3 { // 和| x |<1.5*ln2 
			if !sign {
				hi = x - Ln2Hi
				lo = Ln2Lo
				k = 1
			} else {
				hi = x + Ln2Hi
				lo = -Ln2Lo
				k = -1
			}
		} else {
			if !sign {
				k = int(InvLn2*x + 0.5)
			} else {
				k = int(InvLn2*x - 0.5)
			}
			t := float64(k)
			hi = x - t*Ln2Hi // t*Ln2Hi在这里是精确的
			lo = t * Ln2Lo
		}
		x = hi - lo
		c = (hi - x) - lo
	} else if absx < Tiny { // 当| x |<2**-54时，返回x 
		return x
	} else {
		k = 0
	}

	// x现在在主范围
	hfx := 0.5 * x
	hxs := x * hfx
	r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))))
	t := 3 - r1*hfx
	e := hxs * ((r1 - t) / (6.0 - x*t))
	if k == 0 {
		return x - (x*e - hxs) // c是0 
	}
	e = (x*(e-c) - c)
	e -= hxs
	switch {
	case k == -1:
		return 0.5*(x-e) - 0.5
	case k == 1:
		if x < -0.25 {
			return -2 * (e - (x + 0.5))
		}
		return 1 + 2*(x-e)
		y := 1 - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // 将k添加到y的指数
		return y - 1
	}
	if k < 20 {
		t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k 
		y := t - (e - x)
		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // 将k添加到y的指数
		return y
	}
	t = Float64frombits(uint64(0x3ff-k) << 52) // 2**-k 
	y := x - (e + t)
	y++
	y = Float64frombits(Float64bits(y) + uint64(k)<<52) // 将k添加到y的指数
	return y
}
